pacman::p_load(tidyverse, fishualize, patchwork, ggcorrplot)
path1 <-
'~/Git/henriquelaureano.github.io/THESIS/code/coefs/coefs1to100/'
path2 <-
'~/Git/henriquelaureano.github.io/THESIS/code/coefs/coefs101to200/'
path3 <-
'~/Git/henriquelaureano.github.io/THESIS/code/coefs/coefs201to300/'
fixednames <- c('beta1', 'beta2', 'gama1', 'gama2', 'w1', 'w2')
varnames <- c('logs2_1', 'logs2_2', 'logs2_3', 'logs2_4')
cornames <- c(
'rhoZ12', 'rhoZ13', 'rhoZ14', 'rhoZ23', 'rhoZ24', 'rhoZ34'
)
metnames <- c('conv', 'mll')
vdata <- function(data, path)
{
out <-
tibble::as_tibble(read.table(paste0(path, data)))%>%
dplyr::select(-V1)
type <- substring(data, 6, 7)
switch(type,
v1={ out <-
out%>%'colnames<-'(c(
fixednames,
varnames[1:2], cornames[1], metnames)) },
v2={ out <-
out%>%'colnames<-'(c(
fixednames,
varnames[3:4], cornames[6], metnames)) },
v3={ out <-
out%>%'colnames<-'(c(
fixednames,
varnames, cornames[c(1, 6)], metnames)) },
v4={ out <-
out%>%'colnames<-'(c(
fixednames,
varnames, cornames, metnames)) },
stop('Invalid type'))
return(out)
}
v1data <- v2data <- v3data <- v4data <- vector('list', 18)
label <- numeric(18)
l <- 1
for (i in 1:3)
{
for (j in 1:2)
{
for (k in 1:3)
{
label[l] <- paste0('cs', i, '_', j, '_', k)
v1data[[l]] <-
dplyr::bind_rows(
vdata(paste0('coefsv1_', label[l], '.txt'),
path1),
vdata(paste0('coefsv1_', label[l], '.txt'),
path2),
vdata(paste0('coefsv1_', label[l], '.txt'),
path3)
)
v2data[[l]] <-
dplyr::bind_rows(
vdata(paste0('coefsv2_', label[l], '.txt'),
path1),
vdata(paste0('coefsv2_', label[l], '.txt'),
path2)
)
v3data[[l]] <-
dplyr::bind_rows(
vdata(paste0('coefsv3_', label[l], '.txt'),
path1),
vdata(paste0('coefsv3_', label[l], '.txt'),
path2)
)
v4data[[l]] <-
dplyr::bind_rows(
vdata(paste0('coefsv4_', label[l], '.txt'),
path1),
vdata(paste0('coefsv4_', label[l], '.txt'),
path2)
)
l <- l+1
}
}
}
## str(v1data);length(v1data)
## str(v2data);length(v2data)
## str(v3data);length(v3data)
## str(v4data);length(v4data)
cif1 <- setNames(c( 3, 2.6, 2.5, 4, 5, 10), fixednames)
cif2 <- setNames(c(-2, -1.5, 1, 1.5, 3, 4), fixednames)
logvars <- setNames(log(c(1, 0.6, 0.7, 0.9)), varnames)
rhoZs <- setNames(atanh(c(0.1, -0.5, 0.3, 0.3, -0.4, 0.2)), cornames)
vtrue <- function(data.cif1, data.cif2)
{
out <- rep(c(replicate(data.cif1, n=3, simplify=FALSE),
replicate(data.cif2, n=3, simplify=FALSE)), 3)
return(out)
}
v1cif1 <- c(cif1, logvars[1:2], rhoZs[1])
v1cif2 <- c(cif2, logvars[1:2], rhoZs[1])
v1true <- vtrue(v1cif1, v1cif2)
v2cif1 <- c(cif1, logvars[3:4], rhoZs[6])
v2cif2 <- c(cif2, logvars[3:4], rhoZs[6])
v2true <- vtrue(v2cif1, v2cif2)
v3cif1 <- c(cif1, logvars, rhoZs[c(1, 6)])
v3cif2 <- c(cif2, logvars, rhoZs[c(1, 6)])
v3true <- vtrue(v3cif1, v3cif2)
v4cif1 <- c(cif1, logvars, rhoZs)
v4cif2 <- c(cif2, logvars, rhoZs)
v4true <- vtrue(v4cif1, v4cif2)
cerror <- function(data, true, data_label) {
coefs <- data%>%dplyr::select(-c(conv, mll))
error <- coefs
for (i in seq(ncol(coefs))) { error[i] <- true[i]-coefs[i] }
out <- error%>%
dplyr::summarize_all(mean)%>%
dplyr::add_row(
error%>%
dplyr::summarize_all(quantile, c(0.025, 0.975)))%>%
dplyr::add_row(
error%>%
dplyr::summarize_all(sd))%>%
dplyr::mutate(label=c('mean', 'q025', 'q975', 'sd'),
conv=1-mean(data$conv),
nr=nrow(data),
data=data_label)
return(out)
}
label <- sub('cs1_', 'cs02-', label)
label <- sub('cs2_', 'cs05-', label)
label <- sub('cs3_', 'cs10-', label)
label <- sub('1_', 'high-', label)
label <- sub('2_', '-low-', label)
label <- sub('-1$', '-05k', label)
label <- sub('-2$', '-30k', label)
label <- sub('-3$', '-60k', label)
v1error <- v2error <- v3error <- v4error <- NULL
for (i in seq(18))
{
v1error <-
v1error%>%
dplyr::bind_rows(cerror(v1data[[i]], v1true[[i]], label[i]))
v2error <-
v2error%>%
dplyr::bind_rows(cerror(v2data[[i]], v2true[[i]], label[i]))
v3error <-
v3error%>%
dplyr::bind_rows(cerror(v3data[[i]], v3true[[i]], label[i]))
v4error <-
v4error%>%
dplyr::bind_rows(cerror(v4data[[i]], v4true[[i]], label[i]))
}
## v1error
## v2error
## v3error
## v4error
inside <- function(par, errordata)
{
out <-
errordata%>%
dplyr::select(par, label, conv, nr, data)%>%
tidyr::pivot_wider(names_from=label, values_from=par)
out$data <-
forcats::fct_relevel(forcats::as_factor(out$data), rev)
return(out)
}
biasformat <- function(par, errordatas)
{
out <- dplyr::bind_rows(inside(par, errordatas[[1]]),
inside(par, errordatas[[2]]),
inside(par, errordatas[[3]]),
inside(par, errordatas[[4]]),
.id='v')
out$v <- forcats::fct_recode(out$v,
'RISK MODEL' ='1',
'TIME MODEL' ='2',
'BLOCK-DIAG MODEL'='3',
'COMPLETE MODEL' ='4')
return(out)
}
bias2plot <- function(df, title)
{
ggplot(df, aes(ymin=q025, ymax=q975, x=data, color=conv, size=nr))+
geom_linerange(position=position_dodge(width=0.2))+
geom_point(mapping=aes(x=data, y=mean), size=3,
shape=21, color='black', fill='white', stroke=2)+
facet_grid(~v, scales='free_x')+
geom_hline(yintercept=0, linetype='dashed')+
labs(x=NULL, y='Bias',
title=title,
subtitle='with 2.5% and 97.5% quantiles',
color='Convergance\n(%)',
size='Number of\nmodels')+
coord_flip()+
## theme_minimal()+
## scale_color_viridis_c()+
## scale_color_distiller(palette='Spectral')+
scale_color_fish(option='Trimma_lantana')+
guides(color=guide_colorbar(barwidth=20))+
theme(
strip.background=element_rect(colour='black', fill='white'),
legend.position='bottom',
plot.title=element_text(size=15),
plot.subtitle=element_text(size=14, face='italic'),
axis.text.x=element_text(size=12),
axis.text.y=element_text(size=13),
axis.title.x=element_text(size=13,
margin=unit(c(t=3, r=0, b=0, l=0),
'mm')),
legend.justification=c(1, 0),
legend.title=element_text(size=13),
legend.text=element_text(size=12),
strip.text.x=element_text(size=13, face='bold')
)
}
bias <- function(par, title)
{
out <- biasformat(
par,
list(v1error, v2error, v3error, v4error)
)
bias2plot(out, title=title)
}
bias('beta1', expression(bold('Parameter:'~beta[1])))
bias('beta2', expression(bold('Parameter:'~beta[2])))
bias('gama1', expression(bold('Parameter:'~gamma[1])))
bias('gama2', expression(bold('Parameter:'~gamma[2])))
bias('w1', expression(bold('Parameter:'~w[1])))
bias('w2', expression(bold('Parameter:'~w[2])))
par <- 'logs2_1'
out <- dplyr::bind_rows(inside(par, v1error),
inside(par, v3error),
inside(par, v4error),
.id='v')
out$v <- forcats::fct_recode(out$v,
'RISK MODEL' ='1',
'BLOCK-DIAG MODEL'='2',
'COMPLETE MODEL' ='3')
bias2plot(out, expression(bold('Parameter:'~log(sigma[1]^2))))
par <- 'logs2_2'
out <- dplyr::bind_rows(inside(par, v1error),
inside(par, v3error),
inside(par, v4error),
.id='v')
out$v <- forcats::fct_recode(out$v,
'RISK MODEL' ='1',
'BLOCK-DIAG MODEL'='2',
'COMPLETE MODEL' ='3')
bias2plot(out, expression(bold('Parameter:'~log(sigma[2]^2))))
par <- 'logs2_3'
out <- dplyr::bind_rows(inside(par, v2error),
inside(par, v3error),
inside(par, v4error),
.id='v')
out$v <- forcats::fct_recode(out$v,
'TIME MODEL' ='1',
'BLOCK-DIAG MODEL'='2',
'COMPLETE MODEL' ='3')
bias2plot(out, expression(bold('Parameter:'~log(sigma[3]^2))))
par <- 'logs2_4'
out <- dplyr::bind_rows(inside(par, v2error),
inside(par, v3error),
inside(par, v4error),
.id='v')
out$v <- forcats::fct_recode(out$v,
'TIME MODEL' ='1',
'BLOCK-DIAG MODEL'='2',
'COMPLETE MODEL' ='3')
bias2plot(out, expression(bold('Parameter:'~log(sigma[4]^2))))
par <- 'rhoZ12'
out <- dplyr::bind_rows(inside(par, v1error),
inside(par, v3error),
inside(par, v4error),
.id='v')
out$v <- forcats::fct_recode(out$v,
'RISK MODEL' ='1',
'BLOCK-DIAG MODEL'='2',
'COMPLETE MODEL' ='3')
bias2plot(out, expression(bold('Parameter:'~z(rho[12]))))
par <- 'rhoZ34'
out <- dplyr::bind_rows(inside(par, v2error),
inside(par, v3error),
inside(par, v4error),
.id='v')
out$v <- forcats::fct_recode(out$v,
'TIME MODEL' ='1',
'BLOCK-DIAG MODEL'='2',
'COMPLETE MODEL' ='3')
bias2plot(out, expression(bold('Parameter:'~z(rho[34]))))
bias2plot_rho <- function(df, title)
{
ggplot(df, aes(ymin=q025, ymax=q975, x=data, color=conv, size=nr))+
geom_linerange(position=position_dodge(width=0.2))+
geom_point(mapping=aes(x=data, y=mean), size=3,
shape=21, color='black', fill='white', stroke=2)+
facet_grid(~v, scales='free_x', labeller=label_parsed)+
geom_hline(yintercept=0, linetype='dashed')+
labs(x=NULL, y='Bias',
title=title,
subtitle='with 2.5% and 97.5% quantiles',
color='Convergance\n(%)',
size='Number of\nmodels')+
coord_flip()+
## theme_minimal()+
## scale_color_viridis_c()+
scale_color_fish(option='Trimma_lantana')+
guides(color=guide_colorbar(barwidth=12.5))+
theme(
strip.background=element_rect(colour='black', fill='white'),
legend.position='bottom',
plot.title=element_text(size=15),
plot.subtitle=element_text(size=14, face='italic'),
axis.text.x=element_text(size=12),
axis.text.y=element_text(size=13),
axis.title.x=element_text(size=13,
margin=unit(c(t=3, r=0, b=0, l=0),
'mm')),
legend.justification=c(1, 0),
legend.title=element_text(size=13),
legend.text=element_text(size=12),
strip.text.x=element_text(size=13)
)
}
out <- dplyr::bind_rows(inside('rhoZ13', v4error),
inside('rhoZ24', v4error),
inside('rhoZ14', v4error),
inside('rhoZ23', v4error),
.id='v')
out$v <- forcats::fct_recode(out$v,
'bold(z(rho[13]))'='1',
'bold(z(rho[24]))'='2',
'bold(z(rho[14]))'='3',
'bold(z(rho[23]))'='4')
bias2plot_rho(out,
expression(bold("Complete model's cross-correlations")))
bias2plot2 <- function(df, title)
{
ggplot(df, aes(ymin=mean-1.96*sd, ymax=mean+1.96*sd, x=data,
color=conv, size=nr))+
geom_linerange(position=position_dodge(width=0.2))+
geom_point(mapping=aes(x=data, y=mean), size=3,
shape=21, color='black', fill='white', stroke=2)+
facet_grid(~v, scales='free_x')+
geom_hline(yintercept=0, linetype='dashed')+
labs(x=NULL, y='Bias',
title=title,
subtitle='with \u00b1 1.96 standard deviations',
color='Convergance\n(%)',
size='Number of\nmodels')+
coord_flip()+
## theme_minimal()+
## scale_color_viridis_c()+
## scale_color_distiller(palette='Greys', direction=1)+
scale_color_fish(option='Trimma_lantana')+
guides(color=guide_colorbar(barwidth=20))+
theme(
strip.background=element_rect(colour='black', fill='white'),
legend.position='bottom',
plot.title=element_text(size=15),
plot.subtitle=element_text(size=14, face='italic'),
axis.text.x=element_text(size=12),
axis.text.y=element_text(size=13),
axis.title.x=element_text(size=13,
margin=unit(c(t=3, r=0, b=0, l=0),
'mm')),
legend.justification=c(1, 0),
legend.title=element_text(size=13),
legend.text=element_text(size=12),
strip.text.x=element_text(size=13, face='bold')
)
}
bias2 <- function(par, title)
{
out <- biasformat(
par,
list(v1error, v2error, v3error, v4error)
)
bias2plot2(out, title=title)
}
bias2('beta1', expression(bold('Parameter:'~beta[1])))
bias2('beta2', expression(bold('Parameter:'~beta[2])))
bias2('gama1', expression(bold('Parameter:'~gamma[1])))
bias2('gama2', expression(bold('Parameter:'~gamma[2])))
bias2('w1', expression(bold('Parameter:'~w[1])))
bias2('w2', expression(bold('Parameter:'~w[2])))
par <- 'logs2_1'
out <- dplyr::bind_rows(inside(par, v1error),
inside(par, v3error),
inside(par, v4error),
.id='v')
out$v <- forcats::fct_recode(out$v,
'RISK MODEL' ='1',
'BLOCK-DIAG MODEL'='2',
'COMPLETE MODEL' ='3')
bias2plot2(out, expression(bold('Parameter:'~log(sigma[1]^2))))
par <- 'logs2_2'
out <- dplyr::bind_rows(inside(par, v1error),
inside(par, v3error),
inside(par, v4error),
.id='v')
out$v <- forcats::fct_recode(out$v,
'RISK MODEL' ='1',
'BLOCK-DIAG MODEL'='2',
'COMPLETE MODEL' ='3')
bias2plot2(out, expression(bold('Parameter:'~log(sigma[2]^2))))
par <- 'logs2_3'
out <- dplyr::bind_rows(inside(par, v2error),
inside(par, v3error),
inside(par, v4error),
.id='v')
out$v <- forcats::fct_recode(out$v,
'TIME MODEL' ='1',
'BLOCK-DIAG MODEL'='2',
'COMPLETE MODEL' ='3')
bias2plot2(out, expression(bold('Parameter:'~log(sigma[3]^2))))
par <- 'logs2_4'
out <- dplyr::bind_rows(inside(par, v2error),
inside(par, v3error),
inside(par, v4error),
.id='v')
out$v <- forcats::fct_recode(out$v,
'TIME MODEL' ='1',
'BLOCK-DIAG MODEL'='2',
'COMPLETE MODEL' ='3')
bias2plot2(out, expression(bold('Parameter:'~log(sigma[4]^2))))
par <- 'rhoZ12'
out <- dplyr::bind_rows(inside(par, v1error),
inside(par, v3error),
inside(par, v4error),
.id='v')
out$v <- forcats::fct_recode(out$v,
'RISK MODEL' ='1',
'BLOCK-DIAG MODEL'='2',
'COMPLETE MODEL' ='3')
bias2plot2(out, expression(bold('Parameter:'~z(rho[12]))))
par <- 'rhoZ34'
out <- dplyr::bind_rows(inside(par, v2error),
inside(par, v3error),
inside(par, v4error),
.id='v')
out$v <- forcats::fct_recode(out$v,
'TIME MODEL' ='1',
'BLOCK-DIAG MODEL'='2',
'COMPLETE MODEL' ='3')
bias2plot2(out, expression(bold('Parameter:'~z(rho[34]))))
bias2plot2_rho <- function(df, title)
{
ggplot(df, aes(ymin=mean-1.96*sd, ymax=mean+1.96*sd, x=data,
color=conv, size=nr))+
geom_linerange(position=position_dodge(width=0.2))+
geom_point(mapping=aes(x=data, y=mean), size=3,
shape=21, color='black', fill='white', stroke=2)+
facet_grid(~v, scales='free_x', labeller=label_parsed)+
geom_hline(yintercept=0, linetype='dashed')+
labs(x=NULL, y='Bias',
title=title,
subtitle='with \u00b1 1.96 standard deviations',
color='Convergance\n(%)',
size='Number of\nmodels')+
coord_flip()+
## theme_minimal()+
## scale_color_viridis_c()+
scale_color_fish(option='Trimma_lantana')+
guides(color=guide_colorbar(barwidth=12.5))+
theme(
strip.background=element_rect(colour='black', fill='white'),
legend.position='bottom',
plot.title=element_text(size=15),
plot.subtitle=element_text(size=14, face='italic'),
axis.text.x=element_text(size=12),
axis.text.y=element_text(size=13),
axis.title.x=element_text(size=13,
margin=unit(c(t=3, r=0, b=0, l=0),
'mm')),
legend.justification=c(1, 0),
legend.title=element_text(size=13),
legend.text=element_text(size=12),
strip.text.x=element_text(size=13)
)
}
out <- dplyr::bind_rows(inside('rhoZ13', v4error),
inside('rhoZ24', v4error),
inside('rhoZ14', v4error),
inside('rhoZ23', v4error),
.id='v')
out$v <- forcats::fct_recode(out$v,
'bold(z(rho[13]))'='1',
'bold(z(rho[24]))'='2',
'bold(z(rho[14]))'='3',
'bold(z(rho[23]))'='4')
bias2plot2_rho(out,
expression(bold("Complete model's cross-correlations")))
cmean <- function(data, data_label)
{
df <- data%>%dplyr::select(beta1, beta2, gama1, gama2, w1, w2)
out <-
df%>%
dplyr::summarize_all(mean)%>%
dplyr::mutate(data=data_label)
return(out)
}
lcmean <- function(datas, data_labels)
{
out <- NULL
for (i in seq(18))
{
out <- out%>%bind_rows(cmean(datas[[i]], data_labels[[i]]))
}
return(out)
}
df_cmean <-
dplyr::bind_rows(lcmean(v1data, label),
lcmean(v2data, label),
lcmean(v3data, label),
lcmean(v4data, label), .id='v')%>%
dplyr::mutate(v=v%>%
forcats::fct_recode('RISK MODEL' ='1',
'TIME MODEL' ='2',
'BLOCK-DIAG MODEL'='3',
'COMPLETE MODEL' ='4'))
time <- seq(30, 79.9, length.out=100)
delta <- 80
compCifs <- function(df_cmean.obj)
{
beta1 <- unlist(df_cmean.obj%>%dplyr::select(beta1))
n <- length(beta1)
beta2 <- unlist(df_cmean.obj%>%dplyr::select(beta2))
gama1 <- unlist(df_cmean.obj%>%dplyr::select(gama1))
gama2 <- unlist(df_cmean.obj%>%dplyr::select(gama2))
w1 <- unlist(df_cmean.obj%>%dplyr::select(w1))
w2 <- unlist(df_cmean.obj%>%dplyr::select(w2))
cif1 <- vector(mode='list', length=n)
cif2 <- vector(mode='list', length=n)
for (i in seq(n))
{
risk1 <- exp(beta1[i])
risk2 <- exp(beta2[i])
level <- 1+risk1+risk2
gt <- atanh(2*time/delta-1)
traje1 <- pnorm(w1[i]*gt-gama1[i])
traje2 <- pnorm(w2[i]*gt-gama2[i])
cif1[[i]] <- risk1/level*traje1
cif2[[i]] <- risk2/level*traje2
}
return(list(cif1=cif1, cif2=cif2))
}
cifsHigh <- compCifs(dplyr::filter(df_cmean, grepl('high', data)))
cifsLow <- compCifs(dplyr::filter(df_cmean, grepl( 'low', data)))
getCIF <- function(compCifs.obj, n)
{
cif1 <- cif2 <- NULL
for (i in seq(n))
{
cif1 <- c(cif1, unlist(compCifs.obj[['cif1']][i]))
cif2 <- c(cif2, unlist(compCifs.obj[['cif2']][i]))
}
return(cbind(cif1, cif2))
}
cifs_high <- getCIF(cifsHigh, 36)
cifs_low <- getCIF(cifsLow, 36)
times <- rep(time, 36)
model <- rep(c('RISK MODEL',
'TIME MODEL',
'BLOCK-DIAG MODEL',
'COMPLETE MODEL'), each=900)
labels <- unique(sub('-high|--low', '', label))
labels <- rep(rep(labels, each=100), 4)
dfcif_type <- function(data)
{
out <- tibble::tibble(cif =data,
time =times,
model=forcats::as_factor(model),
label=labels)
return(out)
}
dfcif1_high <- dfcif_type(cifs_high[, 'cif1'])
dfcif2_high <- dfcif_type(cifs_high[, 'cif2'])
dfcif1_low <- dfcif_type(cifs_low[, 'cif1'])
dfcif2_low <- dfcif_type(cifs_low[, 'cif2'])
cif2plot <- function(dfcif_type, descr)
{
ggplot(dfcif_type, aes(x=time, y=cif, color=label))+
geom_line(size=1.5)+
facet_grid(~model)+
theme_minimal()+
labs(x='Time', y=NULL, color=NULL,
title=descr, subtitle='True curve in dashed black')+
## scale_color_fish_d(option='Scarus_quoyi')+
scale_color_brewer(palette='Spectral')+
theme(
axis.text.x=element_text(size=13),
axis.text.y=element_text(size=13),
## legend.position=c(0.05, 0.65),
legend.position='top',
legend.text=element_text(size=14),
legend.box.background=element_rect(color='black'),
strip.background=element_rect(colour='black', fill='white'),
strip.text.x=element_text(size=14, face='bold'),
plot.title=element_text(size=15, face='bold'),
plot.subtitle=element_text(size=15, face='italic'),
axis.title.x=element_text(size=14,
margin=unit(c(t=3, r=0, b=0, l=0),
'mm'))
)+
guides(color=guide_legend(nrow=3))
}
truefixed <-
tibble::tribble(
~beta1, ~beta2, ~gama1, ~gama2, ~w1, ~w2,
3, 2.6, 2.5, 4, 5, 10,
-2, -1.5, 1, 1.5, 3, 4,
)
ciftrue_high <- compCifs(truefixed[1, ])
ciftrue_low <- compCifs(truefixed[2, ])
ciftrue <- function(data)
{
out <- tibble::tibble(time=time, cif=unlist(data))
return(out)
}
cif1true_high <- ciftrue(ciftrue_high[['cif1']])
cif2true_high <- ciftrue(ciftrue_high[['cif2']])
cif1true_low <- ciftrue(ciftrue_low[['cif1']])
cif2true_low <- ciftrue(ciftrue_low[['cif2']])
p1 <-
cif2plot(dfcif1_high, descr='CIF of failure cause 1')+
geom_line(cif1true_high, mapping=aes(x=time, y=cif),
color='black', size=1.5, linetype='dashed')
p2 <-
cif2plot(dfcif2_high, descr='CIF of failure cause 2')+
geom_line(cif2true_high, mapping=aes(x=time, y=cif),
color='black', size=1.5, linetype='dashed')+
theme(legend.position='none')
p1/p2+plot_annotation(title='HIGH CIF SCENARIO',
theme=theme(plot.title=element_text(face='bold')))
p3 <-
cif2plot(dfcif1_low, descr='CIF of failure cause 1')+
geom_line(cif1true_low, mapping=aes(x=time, y=cif),
color='black', size=1.5, linetype='dashed')
p4 <-
cif2plot(dfcif2_low, descr='CIF of failure cause 2')+
geom_line(cif2true_low, mapping=aes(x=time, y=cif),
color='black', size=1.5, linetype='dashed')+
theme(legend.position='none')
p3/p4+plot_annotation(title='LOW CIF SCENARIO',
theme=theme(plot.title=element_text(face='bold')))
COR
coefs22_cor <- round(cor(coefs22), 1)
coefs36_cor <- round(cor(coefs36), 1)
coefs38_cor <- round(cor(coefs38), 1)
coefs40_cor <- round(cor(coefs40), 1)
cor2plot <- function(cor.data, title, zrhos) {
ggcorrplot(cor.data,
type='lower',
lab=TRUE,
## lab_size=5,
## outline.color='white',
colors=RColorBrewer::brewer.pal(n=3, name='Greys'),
title=title,
p.mat=cor_pmat(cor.data),
insig='blank',
ggtheme=theme(
plot.title=element_text(size=16, face='bold'),
axis.text.x=element_text(size=15),
axis.text.y=element_text(size=15),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
panel.background=element_blank(),
axis.ticks=element_blank(),
legend.justification=c(1, 0),
legend.position=c(0.4, 0.85),
legend.direction='horizontal',
legend.title=element_blank(),
legend.text=element_text(size=11)))+
scale_x_discrete(labels=c(expression(beta[2]),
expression(gamma[1]),
expression(gamma[2]),
expression(w[1]),
expression(w[2]),
expression(log(sigma[1]^2)),
expression(log(sigma[2]^2)),
expression(log(sigma[3]^2)),
expression(log(sigma[4]^2)),
zrhos))+
scale_y_discrete(labels=c(expression(beta[1]),
expression(beta[2]),
expression(gamma[1]),
expression(gamma[2]),
expression(w[1]),
expression(w[2]),
expression(log(sigma[1]^2)),
expression(log(sigma[2]^2)),
expression(log(sigma[3]^2)),
expression(log(sigma[4]^2)),
head(zrhos, -1)))+
coord_cartesian()
}
p1 <- cor2plot(coefs22_cor,
title='model22 parameters correlation',
zrhos=c(expression(z(rho[12])),
expression(z(rho[34])),
expression(z(rho[13])),
expression(z(rho[24]))))
p2 <- cor2plot(coefs36_cor,
title='model36 parameters correlation',
zrhos=c(expression(z(rho[12])),
expression(z(rho[34])),
expression(z(rho[13])),
expression(z(rho[24])),
expression(z(rho[14])),
expression(z(rho[23]))))
p3 <- cor2plot(coefs38_cor,
title='model38 parameters correlation',
zrhos=c(expression(z(rho[12])),
expression(z(rho[34])),
expression(z(rho[14])),
expression(z(rho[23]))))
p4 <- cor2plot(coefs40_cor,
title='model40 parameters correlation',
zrhos=c(expression(z(rho[13])),
expression(z(rho[24])),
expression(z(rho[14])),
expression(z(rho[23]))))
(p1|p2)/(p3|p4)
VCOV
plotS <- function(S, title) {
longS <- reshape2::melt(S, na.rm=TRUE)
longS$Var1 <- as_factor(longS$Var1)
longS$Var2 <- fct_relevel(as_factor(longS$Var2), rev)
ggplot(longS, aes(x=Var1, y=Var2, fill=value))+
geom_tile(color='black', size=0.5)+
geom_text(aes(label=round(value, 2)), size=4.5)+
scale_fill_gradient(low='white', high='white')+
labs(title=title)+
scale_x_discrete(labels=c(expression(u[1]),
expression(u[2]),
expression(eta[1]),
expression(eta[2])))+
scale_y_discrete(labels=c(expression(eta[2]),
expression(eta[1]),
expression(u[2]),
expression(u[1])))+
theme_minimal()+
theme(axis.text.x=element_text(size=12, color='black'),
axis.text.y=element_text(size=12, color='black'),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
panel.grid.major=element_blank(),
legend.position='none',
plot.title=element_text(size=12, face='bold'))
}
p1 <- plotS(S=matrix(c(0.4, 0.15, 0.05, 0,
NA, 0.4, 0, 0.05,
NA, NA, 0.25, 0.1,
NA, NA, NA, 0.25), nrow=4, ncol=4),
title='model22: true values')
p2 <- plotS(S=matrix(c(0.2, 0.15, 0.1, 0,
NA, 0.3, 0, 0.1,
NA, NA, 0.4, 0.15,
NA, NA, NA, 0.5), nrow=4, ncol=4),
title='model22: initial guess')
means22 <- colMeans(coefs22)
p3 <- plotS(
S=matrix(c(
exp(means22[7]),
tanh(means22[11])*sqrt(exp(means22[7]))*sqrt(exp(means22[8])),
tanh(means22[13])*sqrt(exp(means22[7]))*sqrt(exp(means22[9])),
0,
NA, exp(means22[8]), 0,
tanh(means22[14])*sqrt(exp(means22[8]))*sqrt(exp(means22[10])),
NA, NA, exp(means22[9]),
tanh(means22[12])*sqrt(exp(means22[9]))*sqrt(exp(means22[10])),
NA, NA, NA, exp(means22[10])),
nrow=4, ncol=4),
title='model22: estimates')
p4 <- plotS(S=matrix(c(0.4, 0.15, 0.05, 0.2,
NA, 0.4, 0.2, 0.05,
NA, NA, 0.25, 0.1,
NA, NA, NA, 0.25), nrow=4, ncol=4),
title='model36: true values')
p5 <- plotS(S=matrix(c(0.2, 0.15, 0.1, 0.2,
NA, 0.3, 0.2, 0.1,
NA, NA, 0.4, 0.15,
NA, NA, NA, 0.5), nrow=4, ncol=4),
title='model36: initial guess')
means36 <- colMeans(coefs36)
p6 <- plotS(
S=matrix(c(
exp(means36[7]),
tanh(means36[11])*sqrt(exp(means36[7]))*sqrt(exp(means36[8])),
tanh(means36[13])*sqrt(exp(means36[7]))*sqrt(exp(means36[9])),
tanh(means36[15])*sqrt(exp(means36[7]))*sqrt(exp(means36[10])),
NA, exp(means36[8]),
tanh(means36[16])*sqrt(exp(means36[8]))*sqrt(exp(means36[9])),
tanh(means36[14])*sqrt(exp(means36[8]))*sqrt(exp(means36[10])),
NA, NA, exp(means36[9]),
tanh(means36[12])*sqrt(exp(means36[9]))*sqrt(exp(means36[10])),
NA, NA, NA, exp(means36[10])), nrow=4, ncol=4),
title='model36: estimates')
p7 <- plotS(S=matrix(c(0.4, 0.15, 0, 0.1,
NA, 0.4, 0.1, 0,
NA, NA, 0.25, 0.1,
NA, NA, NA, 0.25), nrow=4, ncol=4),
title='model38: true values')
p8 <- plotS(S=matrix(c(0.2, 0.2, 0, 0.1,
NA, 0.3, 0.1, 0,
NA, NA, 0.4, 0.15,
NA, NA, NA, 0.5), nrow=4, ncol=4),
title='model38: initial guess')
means38 <- colMeans(coefs38)
p9 <- plotS(
S=matrix(c(
exp(means38[7]),
tanh(means38[11])*sqrt(exp(means38[7]))*sqrt(exp(means38[8])),
0,
tanh(means38[13])*sqrt(exp(means38[7]))*sqrt(exp(means38[10])),
NA, exp(means38[8]),
tanh(means38[14])*sqrt(exp(means38[8]))*sqrt(exp(means38[9])),
0,
NA, NA, exp(means38[9]),
tanh(means38[12])*sqrt(exp(means38[9]))*sqrt(exp(means38[10])),
NA, NA, NA, exp(means38[10])), nrow=4, ncol=4),
title='model38: estimates')
p10 <- plotS(S=matrix(c(0.4, 0, 0.05, 0.2,
NA, 0.4, 0.2, 0.05,
NA, NA, 0.25, 0,
NA, NA, NA, 0.25), nrow=4, ncol=4),
title='model40: true values')
p11 <- plotS(S=matrix(c(0.2, 0, 0.1, 0.2,
NA, 0.3, 0.2, 0.1,
NA, NA, 0.4, 0,
NA, NA, NA, 0.5), nrow=4, ncol=4),
title='model40: initial guess')
means40 <- colMeans(coefs40)
p12 <- plotS(
S=matrix(c(
exp(means40[7]), 0,
tanh(means40[11])*sqrt(exp(means40[7]))*sqrt(exp(means40[9])),
tanh(means40[13])*sqrt(exp(means40[7]))*sqrt(exp(means40[10])),
NA, exp(means40[8]),
tanh(means40[14])*sqrt(exp(means40[8]))*sqrt(exp(means40[9])),
tanh(means40[12])*sqrt(exp(means40[8]))*sqrt(exp(means40[10])),
NA, NA, exp(means40[9]), 0,
NA, NA, NA, exp(means40[10])), nrow=4, ncol=4),
title='model40: estimates')
(p1|p2|p3)/(p4|p5|p6)/(p7|p8|p9)/(p10|p11|p12)
S <- matrix(c(1, 3, 5, 6,
NA, 1, 6, 5,
NA, NA, 2, 4,
NA, NA, NA, 2), nrow=4, ncol=4)
longS <- reshape2::melt(S, na.rm=TRUE)
longS$Var1 <- as_factor(longS$Var1)
longS$Var2 <- fct_relevel(as_factor(longS$Var2), rev)
ggplot(longS, aes(x=Var1, y=Var2, fill=value))+
geom_tile(color='black', size=0.5)+
geom_text(aes( label=c('RISK\nLEVEL', 'RISK\nCORRELATION',
'1st ORDER\nRISK/TIME\nINTERACTION',
'2nd ORDER\nRISK/TIME\nINTERACTION',
'RISK\nLEVEL',
'2nd ORDER\nRISK/TIME\nINTERACTION',
'1st ORDER\nRISK/TIME\nINTERACTION',
'TRAJECTORY\nTIME', 'TIME\nCORRELATION',
'TRAJECTORY\nTIME')), size=5)+
scale_fill_gradient(low='white', high='white')+
scale_x_discrete(labels=c(expression(u[1]),
expression(u[2]),
expression(eta[1]),
expression(eta[2])))+
scale_y_discrete(labels=c(expression(eta[2]),
expression(eta[1]),
expression(u[2]),
expression(u[1])))+
theme_minimal()+
theme(axis.text.x=element_text(size=15, color='black'),
axis.text.y=element_text(size=15, color='black'),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
panel.grid.major=element_blank(),
legend.position='none')